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In light of the rapid recent retreat of Arctic sea Ice, a number of 
studies have discussed the possibility of a critical threshold (or "tip- 
ping point") beyond which the ice-albedo feedback causes the Ice 
cover to melt away in an irreversible process. The focus has typically 
been centered on the annual minimum (September) Ice cover, which 
is often seen as particularly susceptible to destabilization by the ice- 
albedo feedback. Here we examine the central physical processes 
associated with the transition from ice-covered to ice-free Arctic 
Ocean conditions. We show that while the ice-albedo feedback pro- 
motes the existence of multiple Ice cover states, the stabilizing ther- 
modynamic effects of sea Ice mitigate this when the Arctic Ocean 
is ice-covered during a sufficiently large fraction of the year. These 
results suggest that critical threshold behavior is unlikely during the 
approach from current perennial sea ice conditions to seasonally ice- 
free conditions. In a further warmed climate, however, we find that 
a critical threshold associated with the sudden loss of the remaining 
wintertime-only sea ice cover may be likely. 
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The retreat of Arctic sea ice during recent decades (1) 
is believed to be augmented by the difference in albedo 
(i.e., reflectivity) between sea ice and exposed ocean waters 
(2). Because bare or snow-covered sea ice is highly reflective 
to solar radiation, the increasing area of open water that is 
exposed as sea ice recedes leads to an increase in absorbed 
solar radiation, thereby contributing to further ice retreat. A 
number of recent studies have discussed the possibility that 
this positive ice-albedo feedback will cause the rapidly de- 
clining annual minimum (September) sea ice cover to cross a 
critical threshold, after which the sea ice will melt back on an 
irreversible trajectory to a seasonally ice- free state (3-9). 

Heuristically, one might expect in a simple annual mean 
picture of the Arctic Ocean that completely ice-covered and 
ice-free stable states could co-exist under the same climate 
forcing. The ice-free state would remain warm due to the 
absorption of most incident solar radiation, whereas the ice- 
covered state would reflect most solar radiation and remain 
below the freezing temperature. In such a picture, these two 
stable states would be separated by an unstable intermediate 
state in which the Arctic Ocean is partially covered by ice and 
absorbs just enough solar radiation such that it remains at the 
freezing temperature: adding a small amount of additional sea 
ice to this unstable state would lead to less solar absorption, 
cooling, and a further extended sea ice cover. If the back- 
ground climate warmed, the unstable state would require an 
increased ice extent to reflect sufficient solar radiation to re- 
main at the freezing temperature. Beyond a critical threshold, 
the background climate would become so warm that the ice- 
covered state would reach the freezing temperature. At this 
point the stable ice-covered state and unstable intermediate 
state would merge and disappear in a saddle-node bifurcation, 
leaving only the warm ice-free state (10-12). This scenario 
suggests that if an ice-covered Arctic Ocean were warmed be- 
yond the bifurcation point, there would be a rapid transition 
to the ice-free state. It would be an irreversible process in the 
sense that the Arctic Ocean would refreeze only after the cli- 



mate had cooled to a second bifurcation point at which even 
an ice-free Arctic Ocean would become sufficiently cold to 
freeze, representing a significantly colder background climate 
than the original point at which the ice disappeared. Thus 
the ice-albedo feedback could, in principle, cause a hysteresis 
loop in the Arctic climate response to warming. 

Here we investigate the central physical processes under- 
lying the possibility of such a bifurcation threshold in future 
sea ice retreat. We illustrate the discussion with a seasonally 
varying model of the Arctic sea ice-ocean-atmosphere climate 
system. 

Arctic sea ice and climate model 

The theory presented here describes the thermal evolution 
of sea ice, ocean mixed layer, and an energy balance at- 
mosphere that is in steady-state with the underlying sur- 
face forcing, including also representations of dynamic sea ice 
export and diffusive atmospheric meridional heat transport. 
The sea ice thermodynamics in this model is an approxima- 
tion of the full heat conduction equation of Maykut & Un- 
tersteiner (13), which provides the thermodynamic basis for 
most current sea ice models. Ice grows during the winter at 
the base, and when the surface reaches the freezing temper- 
ature in summer, ablation occurs at the surface as well as at 
the base. This model produces an observationally consistent 
simulation of the modern Arctic sea ice seasonal cycle us- 
ing a single one-dimensional nonautonomous ordinary differ- 
ential equation with observationally-based seasonally-varying 
parameters. Here we provide a brief summary of the model 
equations, which are fully derived from basic physical princi- 
ples in the Supporting Information. 

The state variable E represents the energy per unit area 
stored in sea ice as latent heat when the ocean is ice-covered 
or in the ocean mixed layer as sensible heat when the ocean 
is ice-free, 

^_\—Lihi _E < [sea ice] 

~ \cmiHrrdTmi E>0 [oceau] 

where Li is the latent heat of fusion for sea ice, hi is the sea 
ice thickness, Cmi is the mixed layer specific heat capacity, and 
Hmi is the mixed layer depth. The ocean mixed layer tem- 
perature is written in terms of departure from the freezing 
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point, Tml = Tml — Tfr, 



where T,ni is tlie ocean mixed layer 
fr is taken to be 0°C. Tlie time evolution 



temperature and Tj 

of E is proportional to the net energy flux 

'^^ [1 - a{E)] Fs {t) - Fo {t) + AFo 
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dt 



-FT{t)T{t, E)+Fb+ voTli-E). 



[2] 



Here the Stefan-Boltzmann equation for outgoing longwave 
radiation has been linearized in the surface temperature de- 
parture from the freezing point, T{t,E) = T{t,E) — Tfr, as 
Fo{t) -j- FT{t)T{t, E), where the parameters also account for 
the effects of a partially opaque atmosphere and atmospheric 
heat flux convergence which is a function the meridional tem- 
perature gradient. The seasonally varying values of Fo{t) and 
Frit) are derived using an atmospheric model that incorpo- 
rates observations of Arctic cloudiness (f4), surface air tem- 
perature south of the Arctic (f5), and atmospheric transport 
into the Arctic (f6). 

The term AFq represents a specified perturbation to the 
surface heat fiux, which is zero by default but can be increased 
to prescribe a warming in the model. Incident surface short- 
wave radiation Fs{t) and basal heat flux Fb are specified at 
central Arctic values (13). The flnal term in equation [2] 
accounts for an observationally-based constant ice export of 
vo =10% yr~^ (17) when ice is present {E < 0), with the 
ramp function TZ{x) defined to equal x when x > and 
when a; < 0. 

The surface temperature T{t, E) can evolve between three 
different regimes, (i) When ice is present {E < 0) and the sur- 
face temperature is below the freezing point {T{t, E) < 0), it is 
calculated from a balance between the heat flux above the ice 
surface and upward heat flux in the ice, — [1 — 0-{E)] Fs{i) -f 
Fo{t)~-AFo+FT{t)T{t,E) = ~kiT{t,E)/hi = kiLiT{t,E)/E. 
(a) When the surface temperature warms to the freezing point 
(T(t, E) =0), it remains at this point while the ice undergoes 
surface ablation, (in) When the ice ablates entirely, the ocean 
mixed layer is represented as a thermodynamic reservoir us- 
ing T{t,E) = Trni = E/{cmiHmi)- Usiug the ramp function 
as a convenient notation for combining cases (i) and (ii), the 
surface temperature can be expressed as 
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E < [sea ice] 
E > [oceanl 



The ocean is represented as either ice-covered or ice-free at 
any given time. To model the gradual transition between these 
regimes in a partially ice-covered Arctic Ocean, the albedo 
varies between values for ice (a?i) and ocean mixed layer (ami) 
with a characteristic smoothness given by the thickness pa- 
rameter ha, 



a{E) = 



tanh 



Lihn 



[4] 



We also consider a partially linearized version of the model 
in which equation [3] is replaced with 

T{t,E) = —^ [5] 

and there is no ice export {vo = 0). This causes the model 
equations to be linear with the exception of the ice-albedo 
feedback [4]. 

Results 

Seasonal cycle. In a seasonally varying Arctic climate, warm- 
ing might be expected to cause the sea ice to initially melt 
back to the point where the entire Arctic Ocean is ice-free 
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Fig. 1. Sea ice seasonal cycle in a warming climate and solar radiation, (a) Sea- 
sonal cycle of stable solutions of the full nonlinear model are illustrated by plotting 
the model state E (energy per unit area in ocean mixed layer sensible heat or sea 
ice latent heat) versus time of year. Four solutions are plotted, each with different 
levels of surface heating Ai^o- 3 perennial ice state (blue curve, Ai^o — 0), sea- 
sonally ice-free states with most of the year ice-covered (lower red curve, Ai*o = 21 
Wm^"^) or most of the year ice-free (upper red curve, Ai^o = 23 Wm""^), and a 
perennially ice-free state (gray curve, AFq = 19 Wm^'^). As described in equation 
[1], when £/ > 0, it represents the mixed layer temperature of an ice-free ocean 
[E = CjyiiHjj^iTmi)- At E = 0, the ocean mixed layer reaches the freezing 
point {Tml — 0°C), and further cooling will cause ice to grow. When i? < 0, 
it represents the sea ice thici<ness {E = —Lihi); note that ice thickness increases 
downward. Model solutions are drawn with thicker lines when the ocean is ice-covered 
and thinner lines when the ocean is ice-free. Solutions are obtained by integrating 
equations [ 2]-[4] with seasonally varying parameter values given in Table SI {Sup- 
porting Information) until the model has converged on a steady-state seasonal 
cycle. The light gray shaded region to the right represents the first months to be- 
come ice-free in a warming climate (demarcated by zero-crossings of the seasonally 
ice-free solution with Ai^o = 21 Wm^'^), while the light gray shaded region to 
the left represents the last months that are ice-covered in a further warmed climate 
(demarcated by zero-crossings of the seasonally ice-free solution with Ai^o = 23 
Wm^'^). (b) Seasonal cycle of incident solar radiation specified in the model based 
on central Arctic surface observations (13), indicating that the first months to be- 
come ice-free in a warming climate (light gray region to right) and the last months 
to be ice-covered in a further warmed climate (light gray region to left) experience 
similar amounts of solar radiation. Note that the radiation curve is asymmetric due 
to seasonal differences in Arctic cloudiness, but the qualitative results presented here 
do not depend on this asymmetry. 



during part of the year, in contrast to the current perennial 
sea ice cover in the central Arctic. Further warming would 
cause the ice-free period to increase until the Arctic Ocean 
becomes perennially ice-free. We study this scenario theo- 
retically by increasing the imposed surface heat flux AFq in 
equations [2]-[4]. In Fig. a, steady-state seasonal cycle so- 
lutions are plotted in regimes with perennial ice cover (blue 
curve), seasonally ice- free conditions (red curves), and peren- 
nially ice- free conditions (gray curve). 

The annual minimum sea ice area and thickness is com- 
monly referred to as "summer" sea ice and the annual maxi- 
mum is commonly referred to as "winter" sea ice. This nomen- 
clature may carry with it the implication that the ice-albedo 
feedback, which depends on the magnitude of the incident 
solar radiation, would be most prominent during the retreat 
of the summer sea ice cover. Indeed, it is often conjectured 
that a critical threshold for the loss of summer Arctic sea ice 
may be more likely than a threshold for the loss of winter 
ice (8). However, as is illustrated by Fig. b, this terminol- 
ogy can be misleading because the ice cover receives a similar 
amount of incident solar radiation during the period of an- 
nual maximum as at annual minimum. The light gray shaded 
regions in Fig. illustrate the key transition periods in the 
state of the Arctic Ocean during the transition from peren- 
nial ice cover to seasonally ice-free conditions (gray region to 
right) and from seasonally ice-free conditions to perennially 
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Fig. 2. Bifurcation diagram for the partially linearized model, where nonlinear sea 
ice thermodynamic effects have been excluded but the ice-albedo feedback has been 
retained (equations [2]-[3], [5]). For each value of the surface heating AFq, the 
model is integrated until it converges on a steady-state seasonal cycle, and the annual 
maximum (upper curve) and annual minimum (lower curve) values of E are plotted. 
Solutions with perennial sea ice cover are indicated in blue, seasonally ice-free solu- 
tions in red, and perennially ice-free solutions in gray. Dashed lines indicate unstable 
solutions, which have been located by constructing an annual Poincare map and find- 
ing the fixed points (i.e., numerically integrating the model for one year starting from 
an array of initial conditions and identifying the solutions with the same value of E 
at the end of the year as the initial condition). The curves have been smoothed with 
a boxcar filter to suppress a small level of noise associated with numerical integration. 
Note that the lines are slightly curved at the two bifurcation points due to the smooth 
albedo transition associated with ha > 0. The vertical axis is labeled as in Fig. a. 

ice-free conditions (gray region to left). Both of these periods 
experience roughly equivalent amounts of incident solar radi- 
ation (Fig. b), with somewhat more solar radiation occurring 
during the period associated with the loss of winter ice (light 
gray region to left). Hence the ice-albedo feedback should be 
expected to be similarly strong during a transition to peren- 
nially ice-free conditions in a very warm climate (i.e., loss of 
winter ice) as during a more imminent possible warming to 
seasonally ice-free conditions (i.e., loss of summer ice). 

Bifurcation thresholds. We begin the bifurcation analysis us- 
ing the partially linearized version of the model (equations 
[2], [4]-[5]) to focus on the effect of albedo in the absence 
of other nonlinearities. In this representation, the Arctic 
Ocean is viewed as a simple radiating thermal reservoir with 
a temperature-dependent albedo, and the model exhibits a 
linear relaxation to a stable solution in each albedo regime. 
As would be expected by analogy with the discussion above 
of an annual mean Arctic Ocean with a variable sea ice edge. 
Fig. illustrates that when AFq becomes sufficiently large for 
the ocean to remain perennially ice-free with a = Omi, an 
unstable seasonally ice-free solution (red dashed curve) ap- 
pears in a saddle-node bifurcation of cycles (for a discussion of 
the theory of bifurcations in periodic systems, see, e.g., Stro- 
gatz (18)). The unstable solution separates stable solutions 
with perennial ice (blue curve) or perennially ice-free condi- 
tions (gray curve). The perennial ice regime collides with the 
unstable seasonally ice-free state and disappears in a second 
saddle-node bifurcation of cycles at the point where AFq be- 
comes sufficiently large that the ice completely melts at the 
time of annual maximum E in the cold stable state. Due 
to there being significant incident solar radiation during both 
the maximum and minimum periods of the seasonal cycle of 
E (Fig. ), the ice-albedo feedback ensures that all seasonally 
ice- free solutions will be unstable (Fig. ). 

When nonlinear sea ice thermodynamic effects are in- 
cluded (equations [2]-[4]), basal ice formation is controlled 
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Fig. 3. Bifurcation diagram for the full nonlinear model (equations [2]-[4]). 
Axes and colors are as described in the caption of Fig. . The inclusion of nonlinear 
sea ice thermodynamic effects stabilizes the model when sea ice is present during a 
sufficiently large fraction of the year, allowing stable seasonally ice-free solutions (red 
solid curves). Under a modest warming (Ai^o — 15 Wm^'^), modeled sea ice 
thickness varies seasonally between 0.9m and 2.2m. Further warming {/\Fq = 20 
Wm^"^) causes the September ice cover to disappear, and the system undergoes a 
smooth transition to seasonally ice-free conditions. When the model is further warmed 
(Ai^O = 23 Wm^^), a saddle-node bifurcation occurs, and the wintertime sea ice- 
cover abruptly disappears in an irreversible process. While the specific values of Ai^o 
at which the transitions occur are sensitive to parameter choices, the qualitative fea- 
tures of Fig. are highly robust to changes in model parameter values {Supporting 
Information Fig. S4). 

by a diffusive vertical heat flux of kiAT/hi, where AT is the 
difference between surface and basal temperatures and the 
base is assumed to be at the freezing point. This causes thin 
ice to grow significantly faster than thick ice (13). It would 
also cause thin ice to experience greater basal ablation dur- 
ing the summer melt season, but the surface temperature only 
warms until it reaches the freezing point (AT = 0) and surface 
melt begins, making the rate of melt less sensitive to thick- 
ness. These two effects, both nonlinear in E, are expressed in 
equation [3] by the —ki/h = kiLi/E term in the denomina- 
tor and the ramp function TZ{x), respectively. The result is an 
increase in the rate of growth for thin ice which is more stabi- 
lizing for thinner ice, as pointed out (19) and applied (20) in 
previous studies. This is in contrast to the state-independent 
linear mixed layer stabilizing term, —FT{t)E / CmiHmi, which 
applies when E > Q (equations [2] and [3]). 

These nonlinearities allow for the existence of a stable sea- 
sonally ice- free solution (Fig. ). When a sufficiently large 
value of AFq is chosen such that the cold solution becomes 
ice- free during a small part of the year, a slight increase in 
temperature would lead to a longer open-water period and a 
thinner seasonal ice cover. Although the increased period of 
open water promotes warming through the ice-albedo feed- 
back, the thinner ice grows significantly faster because of the 
sea ice thermodynamic effects which are nonlinear in E. Dur- 
ing the ice-covered portion of the year, the stability of the 
solution is controlled by this strong nonlinear stabilizing ef- 
fect, but during the ice-free portion of the year it is replaced 
by the weaker linear mixed layer stabilizing term. This causes 
the stabilizing sea ice thermodynamic effects to dominate the 
destabilizing ice-albedo feedback and allow a stable season- 
ally ice-free solution only when there is ice cover during a 
sufficiently long portion of the year. Nonetheless, the ice- 
albedo feedback causes this regime to warm at an increased 
rate in response to increasing heat flux (compare slopes of 
red and blue curves in Fig. ). As the ice-covered fraction of 
the year decreases in a warming climate, the stabilizing ice 
thermodynamic effects become less pronounced in the full an- 
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nual cycle, and a bifurcation occurs when ice covers the Arctic 
Ocean during a sufficiently small fraction of the year to allow 
the ice-albedo feedback to dominate. Hence when the Arctic 
warms beyond this point, the system supports only an ice- free 
solution (Fig. ). 

Discussion 

Comparison with results of otiier models. The theoretical 
treatment presented here is constructed to facilitate simple 
conceptual interpretation, and to this end many processes 
have been neglected. Factors including possible sea ice-cloud 

feedbacks (21-24), the dcijcndcnce of sea ice surface albedo 
on snow and melt pond coverage (25, 26), ocean heat flux 
convergence feedbacks (6, 27), changes in wind-driven ice dy- 
namics (7), and changes in ice rlieology (28) in a thinning ice 
cover (29) could potentially lead to other bifurcation thresh- 
olds or smooth out the threshold investigated here, akin to the 
smoothing of a first order phase transition due to statistical 
fluctuations (33). We are emboldened in our approach, how- 
ever, because behavior consistent with the mechanism pro- 
posed here can be found in the previously published results 
of models with a broad range of complexities, (i) A "toy 
model" which is forced by a step function seasonal cycle pro- 
duced no stable seasonally ice-free solution in the published 
parameter regime (30), but by a slight adjustment of the tun- 
able model parameters one can find a stable seasonally ice- free 
solution which coexists with a stable perennially ice-free so- 
lution {Supporting Information Fig. S5), consistent with the 
findings presented here, (ii) In a variant of the model used 
in this study that is significantly more complex (representing 
the simultaneous evolution of fractional Arctic sea ice cov- 
erage, moan thickness, and surface temperature, as well as 
ocean mixed layer temperature), increasing the level of green- 
house gas forcing leads to a gradual transition to seasonally 
ice-free solutions followed by a bifurcation threshold during 
the transition to perennially ice- free conditions (31), as in 
Fig. . fiii^ Turning to the most complex current climate mod- 
els, about half the coupled atmosphere-ocean global climate 
models used for the most recent IPCC report (32) jjrcdict 
seasonally ice-free Arctic Ocean conditions by the end of the 
21st century, and none predict perennially ice- free conditions 
by the end of the 2l8t century. However, perennially ice-free 
Arctic Ocean conditions occur in two of the model simula- 
tions after CO2 quadrupling. Neither of the models exhibits 
an abrupt transition when the annual minimum (September) 
ice cover disappears, but after further warming one of the 
models abruptly loses its March ice cover when it becomes 
perennially ice-free (27). The physical mechanism presented 
here may help explain this abrupt loss of simulated March ice 
while the simulated September ice receded gradually. 

Conclusions. Our analysis suggests that a sea ice bifurcation 
threshold (or "tipping point") caused by the ice-albedo feed- 



back is not expected to occur in the transition from current 
perennial sea ice conditions to a seasonally ice-free Arctic 
Ocean, but that a bifurcation threshold associated with the 
sudden loss of the remaining seasonal ice cover may occur 
in response to further heating. These results may be inter- 
preted by viewing the state of the Arctic Ocean as comprising 
a full seasonal cycle, which can include ice-covered periods as 
well as ice-free periods. The ice-albedo feedback promotes the 
existence of multiple states, allowing the possibility of abrupt 
transitions in the sea ice cover as the Arctic is gradually forced 
to warm. Because a similar amount of solar radiation is inci- 
dent at the surface during the first months to become ice-free 
in a warming climate as during the final months to lose their 
ice in a further warmed climate, the ice-albedo feedback is 
similarly strong during both transitions. The asymmetry be- 
tween these two transitions is associated with the fundamental 
nonlinearities of sea ice thermodynamic oft'ccts, which make 
the Arctic climate more stable when sea ice is present than 
when the open ocean is exposed. Hence when sea ice covers 
the Arctic Ocean during fewer months of the year, the state of 
the Arctic becomes less stable and more susceptible to desta- 
bilization by the ice-albedo feedback. In a warming climate, 
as discussed above, this causes irreversible threshold behavior 
dming the jjotential distant loss of winter ice, but not during 
the more inunincnt possible loss of summer (September) ice. 

The relevance of any basic theory to the actual future 
evolution of the complex climate system must be carefully 
qualified. Since the time scale associated with the sea ice re- 
sponse to a change in forcing may be decadal, and the time 
scale associated with increasing greenhouse gas concentrations 
may be similar, the system may not be operating close to a 
steady-state. In the gradual approach to steady-state under a 
continual change in forcing, the difi'erence between a region of 
the steady-state solution with increased sensitivity to the forc- 
ing and an actual discontinuous bifurcation threshold (as in 
Fig. ) could be difficult to discern. If greenhouse gas concen- 
trations were reduced after crossing a bifurcation threshold, 
however, the possible irreversibility of the trajectory would 
certainly be expected to be relevant. 
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Here we derive the idealized Arctic sea ice-ocean- 
atmosphere model that is summarized in equations [l]-[4] 
of the Research Report. Note that we carry out the entire 
derivation using dimensional variables, rather than following 
the conventional mathematical development of such equations 
through a process of non-dimensionalization, in order to make 
direct contact with previous studies of the thermodynamics of 
sea ice and climate. 



MU71, whereas in the full coupled version of the model Fo{t) 
and Frit) take on a different set of values computed using the 
atmospheric model (equations [39]-[40] below). 

At the ice-ocean interface (z = hs), MU71 apply a Stefan 
condition for ice growth or ablation, 



dJiB 
'~dr 



kea{T, S) 



dT 

dz 



Fb, 



Sea Ice 

The evolution of the sea ice temperature profile is an ide- 
alized version of the single-column thermodynamic model of 
Maykut and Untersteiner (1) (hereafter MU71). Vertical heat 
conduction in sea ice is computed in MU71 according to 



Ceff(T,S)- = - 



fceff (T, S) 



dT 

dz 



Ai 



which can be derived from the general theory of mushy layers 
(2). Here An represents the absorption of shortwave radiation 
that has penetrated below the surface of the ice, the effective 
heat capacity Cch{T,S) and thermal conductivity keff{T,S) 
depend on simulated temperature T and specified salinity S, 
and the vertical coordinate z increases upward. Note that for 
the T and S range in perennial ice, MU71 neglect the vertical 
derivative of the effective conductivity, dki.g(T, S)/dz, allow- 
ing the first term on the right-hand side of equation [6] to 
be expressed as fceft(T, S)d'^T/dz^. MU71 also include a layer 
of snow above the ice with specified snowfall and simulated 
snow melt. 

The boundary condition in MU71 at the upper surface 
(« = /it) is a fiux balance when the ice is below the freez- 
ing temperature (T/r) and otherwise a Stefan condition for 
surface ablation: 



fceft(T, S) 



dT 

dz 



Ftop{t,Ti,ai) = 



T<0 
Ti = ' 



[r] 

with Li the latent heat of fusion of ice, the surface albedo, 
Ti = Ti — Tfr the surface temperature departure from the 
freezing point with T = T{z = Ht), and Ftop{t,Ti, at) repre- 
senting the sum of sensible and latent heat fluxes and long- 
wave and shortwave radiative fluxes out of the surface. The 
seasonal cycle of each of these components of the surface flux 
are specified in MU71 based on observations, except for the 
upward longwave fiux which is computed from the surface 
temperature using the Stefan-Boltzmann equation. To facil- 
itate an analytical solution for T (equation [15] below), we 
approximate the Stefan-Boltzmann equation by its linearized 
version, crT/ = (Jq + o-tT, where a is the Stefan-Boltzmann 
constant and the parameters [ag = 316 Wm^^, = 3.9 
are chosen such that the equation is exact when 
Ti = — 30°C and when T = 0°C, which are the approximate 
values of T during most of the winter and summer, respec- 
tively. This allows the temperature dependence of the surface 
fiux to be expressed as 

Ftopit,Ti,a,) = -{1 - ai)Fs{t) + Fo{t) + FT{t)T, [8] 

where Fs{t) is the downwelling shortwave radiation fiux, Fo{t) 
is ctq plus the specified sensible and latent heat fluxes, and 
Frit) = fTT. Note that here the atmosphere is specified as in 



with the fiux from the ocean into the base of the ice specified 
to take a constant value of Fb = 2 Wm^^. Note that the 
temperature at the ice-ocean interface must be at the freez- 
ing point, T{z = Hb) = Tfr- The upper and lower surfaces 
of the ice, hr and Hb, evolve separately in MU71, who use a 
coordinate system in which each ice parcel remains stationary, 
and the predicted ice thickness is hi = hr — fiB (see schematic 
in Fig. SI). 

Here we neglect snow (MU71 report that having no snow 
causes the annual mean thickness to increase by 17cm from 
the standard case value of 288cm), and we neglect penetrating 
shortwave radiation Ar (which MU71 report causes the 
annual mean ice thickness to decrease by 45cm). The impact 
of neglecting both of these factors is shown in Fig. S2 (black 
curves). 

The thermal conductivity fcefi(T, S) in the MU71 stan- 
dard case run is always 90%-100% of the pure ice value, 
and we approximate it to take the constant pure ice value, 
fccff (T, 5) = fci = 2 Wm^^K^^. The freezing temperature in 
MU71 is taken to be Tfr = — 1.8°C at the base of the ice and 
Tfr = — 0.1°C at the upper ice surface, and we approximate 
it to take a constant value of Tfr = 0°C. Lastly, because the 
Stefan number Ns = Li/ |^Ceff(T, S')AT'j is large, the temper- 
ature field in the ice relaxes quickly in response to changes 
in the solidification rate. The actual values of Ns predicted 
by the MU71 standard case seasonal cycle vary widely but 
typically are Ns 3> 1. Under these conditions, the system 
[6]-[7], [9] can be expressed as a single ordinary differential 
equation, as described below. 

Applying the large Stefan number approximation to the 
heat conduction equation [6] yields a linear temperature pro- 
file, f = ffr + [T - ffr){z - hB)/{hT - Hb) = Tfr + Ti{z - 

hB)/hi, which can be derived using a scaling argument after 
vertically integrating equation [6] and inserting the boundary 
conditions [7], [9]. Next, this quasi-stationary temperature 
field is inserted into equations [7], [9]. The upper bound- 
ary condition [7] includes two different cases depending on 
whether or not surface ablation is occurring: 



atm 



Z=/7, 



top 



ice 



T=r. 



h. 



ocean 



T=T 



fr 



B 



Fig. SI. Schematic showing fluxes and variables in the sea ice component of 
the idealized model presented here. All fluxes are defined such that a positive value 
implies an upward flux. 
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Fig. S2. Effects of approximating tlie ice thermodynamics in tlie model of MU71 (1). (a) Steady-state solution seasonal cycle of ice thickness in the MU71 standard case 
simulation (gray curve and circles), in a simulation with the MU71 model carried out for this study with no snow or penetrating shortwave radiation {hs = A.ji = 0; black 
curve and circles), when the MU71 representation is replaced by the idealized sea ice model given by equations [15]-[16] (blue curve), and the standard case run with the 
fully coupled idealized sea ice-ocean-atmosphere model summarized in equations [l]-[4] of the Research Report (red curve), (b) Seasonal cycle of surface temperature for 
the same four solutions as in (a). Note that the surface temperature in the idealized model is diagnosed from the computed thickness and the specified surface forcing, (c) 
Relaxation time to reach steady-state ice thickness from two different initial conditions in the MU71 model with hg = Aji = (black curves and circles) and in the idealized 
sea ice model (blue curves). 



(i) When the surface is below the freezing temperature 
(Ti < 0), the upper boundary condition [7] with the hnear 
temperature profile takes the form 



h-r- = —Ftop{t,Ti,ai). 
hi 



[101 



Since there is no surface ablation [dhr/dt = 0), the ice thick- 
ness evolves as dhi/dt = d/dt(hT — hs) = —dhs/dt, and af- 
ter inserting the linear temperature profile the lower boundary 
condition [9] becomes 



Li 



dhi 



hi 



hi 



[11] 



Inserting the surface temperature [10] into equation [11] 
shows that ice thickness evolves according to 



dhi 
"'dt 



[12] 



(ii) During surface ablation [Ti = 0), 
profile takes on a constant value of T = Tfr. 
boundary condition [7] takes the form 

dhr 

Ftop{t,Ti,ai) = Li-^, 
and the lower boundary condition [9] becomes 

dt 

which together imply that ice thickness, hi = hr — hB, evolves 
according to an equation identical to the case with the surface 
below the freezing temperature [12]. 

The steady-state surface temperature can be derived by 
inserting equation [8] into equation [10] to yield an algebraic 
solution for the case with Ti < 0, which can be combined with 
the ablation case {Ti = 0) as 



the temperature 
Hence the upper 

[13] 



[141 



Ti{t,hi) 



{l-ai)Fs{t)- Fo{t) 
-ki/hi-Frit) 



[15] 



Here the dependence of T on t and hi has been explicitly indi- 
cated, and the ramp function TZ{x) is defined to be TZ{x) = 
if X < and 7?.(a::) = a::ifa;>0. Note that the two surface 
boundary conditions in equation [7] are compactly embodied 
in the ramp function in equation [15]. The thickness evolu- 
tion in both cases [12] can be written after inserting equation 
[8] as 

Li^ = -{l-ai)Fs{t) + FQ{t)+FT{t)Ti{t,hi)-FB. [16] 
at 



The sea ice model is fully contained in equations [15]- 
[16]. The results of this idealized ice thermodynamics model 
forced by specified surface and basal fluxes as in MU71 are 
shown in Fig. S2 (blue curves), which indicates that this 
approximate representation yields results in good agreement 
with the full numerical solution to the partial differential equa- 
tion [6] in MU71 (cf. refs. 3, 4). 

While most aspects of horizontal sea ice dynamics are ne- 
glected in this idealized treatment, in the coupled version of 
the model (equations [l]-[4] of the Research Report) we pa- 
rameterize the net annual export of sea ice out of the central 
Arctic, most of which escapes through Fram Straight. Arctic 
sea ice has a residence time of roughly 3-12 years (5), with 
a net annual export of about 10% of the ice area (6). This 
continuous export makes the ice thickness somewhat more sta- 
ble: to maintain thicker ice, a larger amount of new ice must 
be produced each year. We approximately account for this 
by adding to the ice thickness evolution [16] a decay term 
—voLihi, with vo = 0.1 yr^^. 



Atmosphere 

In the presence of significantly different Arctic Ocean surface 
conditions, such as an exposed ocean mixed layer, the atmo- 
spheric energy fluxes into the surface are also expected to 
change significantly. This is particularly true for the down- 
welling longwave radiation which includes the effects of both 
horizontal atmospheric heat flux convergence and downward 
emission of absorbed upward longwave radiation due to the 
opacity of the atmosphere (i.e., the greenhouse effect). Here 
we use an idealized atmospheric model to account for changes 
in downwelling longwave radiation. This allows us to ap- 
proximate Ftop{t,T,a) over a wide range of climates. The 
derivation that follows is similar to previous treatments of 
two-stream radiative atmospheres (e.g., refs. 3,7,8). 

Heat Flux Convergence. The meridional heat flux convergence 
averaged over 70°N-90°N is equivalent to a spatially averaged 
vertical fiux of roue hly D = 100 Wm"^ (9). Since the pole- 
ward heat flux in the atmosphere is related to transport of 
sensible and latent heat by eddies, it is often approximated in 
idealized climate models as being proportional to the merid- 
ional temperature difference (10, 11, 12), which is equivalent 
to assuming meridional effective diffusion of temperature as 
in typical atmospheric energy balance models (13, 14, 15). Al- 
though a destabilizing increase in atmospheric meridional heat 
flux into the Arctic may occur in response to warming due to 
factors including increased humidity (16, 17, 18, 19), if the 
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warming is significant then reduced atmospheric heat trans- 
port is expected to be a principal damping mechanism (20). 
Here we follow the convention of setting the meridional heat 
flux to be proportional to the meridional temperature differ- 
ence, 

D{t,T) = koATmeridit), [17] 

where ATmerid = Tsouth{t) — T with T the simulated surface 
temperature in the Arctic and Tsouth it) the seasonally varying 
temperature south of the Arctic which is specified here from 
NCEP-NCAR reanalysis 1971-2000 climatological lOOOmb at- 
mospheric temperature (21) spatially averaged from the equa- 
tor to 70°N. We use ko = 2.7 Wm^^/K, which optimizes the 
match to observed poleward heat transport (9) (although this 
parameterization leads to a model annual cycle in D that is 
somewhat exaggerated compared to observations). 

Longwave Absorption. We use a vertically continuous dry en- 
ergy balance atmospheric model. We approximate there to 
be no absorption of shortwave radiation in the atmosphere 
and no scattering of longwave radiation. Longwave radiation 
is absorbed and emitted in continuous vertical levels with an 
absorption cross section that is independent of wavelength, 
temperature, and pressure, and the radiative fields are solved 
using a two-stream approximation. We assume that the pole- 
ward atmospheric heat transport into the Arctic, D{t,T), is 
distributed uniformly in optical height (3). 

The intensity of a beam of radiation propagating verti- 
cally upward from Earth's surface will diminish with height 
z according to dl/dz = —p{z)k{z)I, where p{z) is the atmo- 
spheric density and k{z) is the extinction coefficient. This can 
be solved for intensity as a function of height. 



I = lo exp 



p{z')K{z')dz' 



= lo exp[— r(2) 



[181 



where Iq is the intensity at the surface and the optical height 
is r{z) = p{z')K{z')dz' . We measure height using t{z) in- 
stead of z, which has the advantage that k{z) and p{z) are 
eliminated from the equations and the atmosphere can be ap- 
proximately described by a single parameter, the total opti- 
cal thickness Ti = t{z — > oo). Note that our use of optical 
height differs slightly from the standard convention of using 
optical depth, which is integrated from the top of the atmo- 
sphere downward. Regarding the physical meaning of the op- 
tical thickness ti , note that the fraction of longwave radiation 
emitted vertically from the surface that escapes to space is 
exp(— Ti). A slanted optical path in the atmosphere, St* , can 
be related to optical height according to 6t = 6t* cos 6, where 
9 is the angle the path makes with the vertical (Fig. S3). 

The intensity of longwave radiation in the atmosphere, 
I{T,d,(j)), is a scalar field that depends on optical height and 
direction, with 4> being the azimuthal angle (Fig. S3). We 
model the atmosphere as a grey material that absorbs a frac- 
tion St* of the intensity passing through it and emits an equiv- 
alent fraction St* of its blackbody radiation. The blackbody 
radiation of an air parcel can be computed from the Stefan- 
Boltzmann equation, B(Ta) = aT^/ir, where Ta is the temper- 
ature of the air parcel and the factor n accounts for radiation 
occurring in every direction from a point source. This gives 
a change in intensity of 51 = —I St* + BSt* , which becomes 
the Schwarzschild equation when written in terms of optical 
height: 

^^^("'^''^)=-7(r,e,0) + B(r). 



cos 6- 



[191 



Note that B is independent of angle since a blackbody emits 
radiation equally in all directions. 



We assume horizontally uniform radiation from the sur- 
face and a horizontally homogenous atmospheric medium, 
which makes the intensity horizontally isotropic due to az- 
imuthal symmetry, 7(r, 6, 0) = I(t,9). In thermal steady- 
state, this can be written as a vertically constant divergence 
of vertical net flux. 



d_ 

dT 



dcj7(r, 9) cos 9 = 



D{t,T) 



[201 



Ti 



where we have deflned the integral over all solid angles, 
f dio = d9 f^^ sin 9d4>. We can rewrite the divergence con- 
dition [20] as an algebraic equation by equating it with the in- 
tegral over all solid angles of the Schwarzschild equation [19]. 
Solving this for B{t) and inserting this into the Schwarzschild 
equation [19] gives a single integro-differential equation for 
I{r,9), 



cos 8 



dI{T, 

dT 



= -I{t,9) 



1 

47r 



D{t,T) 



dujI{T, 9) 



21] 

The boundary conditions are that there is no downward 
longwave radiation at the top of the atmosphere (r = ri). 



7(ri,7r/2 < 9 < n) = 0, 



[221 



and that the upward radiative flux from the surface (r = 0) 
is (70 -I- arT, 



7(0,0 < 9 < 71-/2) = 



o^o + (JtT 



[231 



The system [21]- [23] 
wave radiation 7(r, 9) 



uniquely specifies atmospheric long- 
given the poleward heat transport 
D(t,T), atmospheric optical thickness ri, and surface tem- 
perature T. 

Here we use a standard two-stream approximation to ar- 
rive at an analytical solution to the system [21]-[23]. Multi- 
plying the Schwarzschild equation [19] by cos 6 and integrat- 
ing over the upper and lower hemispheres (i.e., the upper and 
lower halves of an infinitesimal sphere surrounding a given 
point in the atmosphere) leads to 



d_ 

dT 

d_ 

dT 



du^ I{t, 9) cose = -F\t)+-kB{t), 



[241 



j du^I[T,9)cos9 = F\t)-ttB{t), [25] 

where we have defined total upward and downward fiuxes 
through a horizontal surface in the atmosphere. 



duj^ I {t, 9) COS 9, 



[261 
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Fig. S3. Schematic of atmospheric model for computing Ftop{t,T^ ex.). 
D{t,T) represents meridional heat transport, and (1 — a)Fs amount of 

absorbed solar radiation. Longwave radiative intensity is represented by / and optical 
height is given as r with 5t* an optical path at angle 6 to the vertical. Note that here 
the optical height increases upward, in contrast with the optical depth which increases 
downward and is typically used in radiative transfer calculations. The total upward 
and downward longwave radiative fluxes through a horizontal surface are F-i^{t) and 
F— (r) , respectively. The model allows the surface incident longwave radiation to be 
represented as a function of outgoing surface longwave radiation (equation [34]). 
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duj^ I {r,e) cose, 



[27] 



and the integrals over solid angles in each hemisphere 
are defined as J du)^ = J^^^dO J^^ sin 9 d<f) and J duj^ = 
Xr/2 de f^'' sin ed(f>. 

The type of two-stream approximation we employ, an ex- 
ponential kernel approximation (sec. 2.4 of ref. 7), is equiva- 
lent to assuming that radiation propagates through the atmo- 
sphere only at one effective angle. This angle deviates from 
the vertical to account for the fact that in the exact solution 
radiation propagates in all directions. Upwelling radiation 
propagates along 6 = and downwelling radiation propa- 
gates along 6 = ^off + TJ". This assumption allows equations 
[24]-[25] to be written as 



cos ^eff 



cos 6ea 



dFHj) 
dr 



-F\t)+7tB{t), 
F-^{t) - ttB{t). 



[28] 



[29] 



Note that the negative sign from the definition of [27] is 
cancelled in equation[29] because cos {6cB + vr) = — cos 9cS. 

Inserting the definitions [26]-[27] into the divergence 
condition [20] leads to 



_d_ 
dr 



-F 



D{t,T) 



[30] 



Ti 



which can be multiplied by cos O^a and then equated with 
the sum of equations [28] and [29]. This allows us to write 
the divergence condition as an algebraic relation that can be 
directly solved for B{t): 



B{r) 



1 

2^ 



D{t,T) cos e^ft 



+ F' 



Ti 



[31] 



Finally, inserting the definitions [26]-[27] into the bound- 
ary conditions [22]- [23] leads to 



^1 



0, 



F\0) 



Co + (7tT. 



[32] 



[33] 



Equations [28]-[29], [31] represent a system of two cou- 
pled first order linear ordinary differential equations, which 
can be solved subject to boundary conditions [32]- [33] to 
give the full vertical profiles of F^ (r) and F^ (r) . Note that the 
temperature profile can be calculated from this solution using 
equation [31] and the Stefan-Boltzmann equation. From the 
full solution to the system [28]-[29], [31]-[33] (not shown), 
the surface downward longwave radiation is 



F^{0) 



(TO + (JtT D{t,T) 



[34] 



The incident surface radiation includes half of the atmospheric 
heat convergence, D{t, T)/2, while the other half is emitted to 
space. This arises from the assumption that D{t,T) is evenly 
distributed in optical height. 

The choice of 6eff can be optimized to match the ex- 
act solution of the system [21]-[23] given values of D{t,T), 
Ti, and FT(0) = ffo + (ttT. Salby (sec. 8.4.2 of ref. 8) 
finds that the effect of averaging over spectral bands sug- 
gests the value cos^eff = |. Goody and Yung (sec. 9.2.1 
of ref. 7) derive an atmosphere similar to the system [28]- 
[29], [31], with D{t,T) = 0, by using a version of the two- 
stream approximation in which hemispheric isotropy is as- 
sumed: I{T,e > 7r/2) = 7+ and 7(r, 61 > tt/2) = This 
yields a result equivalent to letting cos^eff = |. Thorndike 



(3) derives an atmosphere analogous to the model derived 
here but with only vertically propagating radiation, arriving 
at a result equivalent to cos^eff = 1. We take the optical 
thickness ri as the tunable parameter in the model, and we 
leave cos^eff unspecified because it simply scales the optical 
thickness. 

The net longwave radiation from the solution [34] is 
F\0)-F^{0) = f,^„{ao+aTT)-^^i^, [35] 
where we have defined an atmospheric greenhouse factor as 



1 



f^LW 



1 + 



1 + 



[36] 



Note that < < 1. This makes clear the effect of the 

atmospheric longwave radiation model: it mitigates surface 
longwave cooling by emitting some of the energy back to the 
surface, equivalent to reducing the surface upward longwave 
radiation by the factor Klw, and it adds energy associated 
with atmospheric heat flux convergence (i.e., net meridional 
heat transport into the Arctic). By reducing the surface net 
upward longwave radiation, the interactive atmospheric model 
weakens the stabilizing influence of outgoing longwave radia- 
tion on the ice/ocean system. 

Grey two-stream atmospheres like the model used here can 
capture many of the basic features of radiative transfer in an 
approximate way. A thorough comparison of various types of 
two-stream approximations is discussed in Goody and Yung 
(sec. 2.4 of ref. 7). 

The optical thickness of the atmosphere depends on water 
vapor, cloud particles, and greenhouse gases such as carbon 
dioxide. It is higher during summer than during winter be- 
cause of increased water vapor and cloudiness. We specify 
the optical thickness seasonal cycle to follow observed Arctic 
cloudiness, 

TO + Tc,fc{t) 

where /c(i) is the Arctic cloud fraction seasonal cycle speci- 
fied from observations (22) and tq and Tc are chosen to give a 
sea ice seasonal cycle matching that computed using forcing 
from MU71 (cf. refs. 3, 23). This leads to a choice of tq = 0.5 
and Tc = 3.6. 

The actual energy fiux at the top of the sea ice or ex- 
posed ocean mixed layer, Ftop{t,T,a), includes components 
of sensible and latent heat fluxes in addition to downward 
and upward shortwave and longwave radiation. According 
to the observationally-based central Arctic values specifled in 
MU71, the sensible and latent heat fluxes are small compared 
to the radiative fluxes, and here we effectively approximate 
the sensible and latent heat fluxes by incorporating them into 
the computed downwelling longwave flux. The longwave emis- 
sivities of ice and open water, both roughly 0.95-1, are here 
approximated to unity. Under these approximations, the total 
surface flux can be written 

Ftop{t,T,a) = Ki„,.(t) {go + otT) - i^Ml - (1 - a)Fs{t), 

[38] 

where the downwelling shortwave radiation incident at the 
surface, Fsit), is specified from observations as in MU71. In- 
serting equation [17], we see that the full temperature de- 
pendence is linear, allowing us to write Ftop [38] in the form 
of equation [8] with parameters 



Fo{t) = Klw{1:)(to ^T'sotit/i(i) 



and 



Frit) = KL„{t)aT + 



ki 



[391 



[401 
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Ocean Mixed Layer 

To allow the simulation of ice-free conditions, we include a rep- 
resentation of an ocean mixed layer which becomes exposed 
when all of the ice ablates. The mixed layer is represented 
as a thermodynamic reservoir with a characteristic depth of 
Hmi = 50m, in agreement with observations (24). The mixed 
layer temperature evolution is proportional to the net flux, 

CmlHral-^ = (1 - QmO Fg (t) - Fq (t) - FT{t)Tml + Fb , 

[41] 

with mixed layer heat capacity Cmi = 4 x lO'' Jm '^K ^. We 
use an open water albedo of ami = 0.2, similar to previous 
studies (3, 25), to account for the presence of small amounts 
of thin ice in a largely ice-free Arctic Ocean. When the ice 
completely melts {hi = 0), the ocean mixed layer temperature 
is evolved forward from Tmi = 0, and when the mixed layer 
cools back to Tmi = 0, the ice thickness is evolved once again 
starting from hi = 0. 



Coupled Model 

The separate equations for Tmi and hi can be combined, since 
only one is evolving at any given time. We define the energy 
per unit area in the system, E, to be equal to the sum of the la- 
tent heat content of the sea ice and the specific heat content of 
the ocean mixed layer (equation [1] in the Research Report). 
This allows the ice and ocean mixed layer components of the 
idealized model [15]- [16], [41] to be expressed as equations 
[2]-[3] in the Research Report. The parameters Fo{t) and 
Frit), which are used to determine the surface energy flux, 
have values computed using the atmospheric model [39]-[40]. 
An imposed annually constant surface energy flux is included 
in equations [2]- [3] by replacing Fo{t) with Fo(t) — AFq. Val- 
ues for the parameters in equations [l]-[4] in the Research 
Report are given in Table SI. The standard case simulation 
with this model, illustrated in Fig. S2, produces central Arctic 
sea ice conditions in fairly good agreement with MU71. 



Table SI. Descriptions and default values of model parameters. Time evolution t is measured in years while fluxes are measured in Wm~^, 
which allows most dimensional parameters to be approximately of order unity but requires a non-conventional choice of units for energy per unit 
area E (written in Wm^^yr), heat capacity CmiHmi, and latent heat Li. For the three seasonally varying parameters, the annual mean value is 
given in the table; the monthly values starting with January are Fo{t) = (120, 120, 130, 94, 64, 61, 57, 54, 56, 64, 82, 110) Wm"^, Frit) = 
(3.1, 3.2, 3.3, 2.9, 2.6, 2.6, 2.6, 2.5, 2.5, 2.6, 2.7, 3.1) Wm-^K-i, and Fs(t) = (0, 0, 30, 160, 280, 310, 220, 140, 59, 6.4, 0, 0) Wm^^. 



Symbol 


Description 


Value 




Latent heat of fusion of ice 


9.5 Wm^^yr 




Ocean mixed layer heat capacity times depth 


6.3 \Nm-^yrK-^ 


Cti 


Albedo when surface is ice-covered 


0.68 


Ctml 


Albedo when ocean mixed layer is exposed 


0.2 


ki 


Ice thermal conductivity 


2 Wm^^K-i 


Fb 


Heat flux into bottom of sea ice or ocean mixed layer 


2 Wm-2 


ha 


Ice thickness range for smooth transition from Oi to ami 


0.5 m 


vo 


Dynamic export of ice from model domain 


0.1 yr-^ 


Fo{t) 


Temperature-independent surface flux (seasonally varying) 


85 Wm"^ 


Frit) 


Temperature-dependent surface flux (seasonally varying) 


2.8 Wm^^K-i 


Fs{t) 


Incident shortwave radiation flux (seasonally varying) 


100 Wm-=^ 


AFo 


Imposed surface heat flux 


Wm-2 
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Fig. S4. Robustness of the results in Fig. 3 of the Research Report to parameter 
regime, illustrated here by varying the parameter governing the smoothness of the 
albedo transition (/ia). For each value of ha_, ranges of surface heating (AFq) 
that give rise to stable solutions that are perennially ice-covered (blue region), sea- 
sonally ice-free (red region), or perennially ice-free (white region) are identified. The 
default parameter regime is indicated by a black diamond. Mixed shades indicate 
the overlap in regions where multiple stable solutions coexist, and the bifurcation 
curves marking the edges of this space are indicated by black curves. The lack of any 
purple region, which would indicate an overlap between red (seasonally ice-free) and 
blue (perennial ice), demonstrates that multiple states are not found with the warm 
state being seasonally ice-free, while the presence of light red and light blue regions 
shows that multiple states with the warm state being perennially ice-free are possible. 
The variation of other model parameters (not shown) leads to similar results. This 
indicates that although the size of the A_Fo range where multiple solutions coexist 
depends on ha_, both (i) the lack of a bifurcation threshold during the transition 
from perennial ice to seasonally ice-free conditions and (ii) the presence of multiple 
states and threshold behavior during the transition to perennially ice-free conditions 
are robust features of the model equations. 
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Fig. S5. Seasonal cycle of Arctic sea ice and ocean conditions simulated with 
the toy model of Thorndike (3). Seasonally varying solutions are plotted as closed 
curves in the three-dimensional model state space, which represents changes in sea 
ice thickness, sea ice surface temperature, and ocean mixed layer temperature. The 
standard case solution is indicated by the black curve. Thorndike found that when the 
specified atmospheric Arctic heat flux convergence D was increased, the model tran- 
sitioned from perennially ice-covered to perennially ice-free conditions, with no stable 
seasonally ice-free solution possible. Rather than prescribing observed seasonally- 
varying forcing quantities, Thorndike assumed a step-function form for the forcing, 
with shortwave radiation and optical thickness taking on constant values during the 
summer and winter half-years. He found observationally consistent ice thickness with 
summer and winter optical thicknesses of 4.5 and 3, respectively (black curve). When 
we choose for these parameters instead 1.5 and 5, respectively, Thorndike's toy model 
simulates a relatively consistent approximation of the modern sea ice seasonal cycle 
(blue curve). Increasing D from 100 Wm""^ to 145 Wm~^ in this regime, however, 
produces a stable seasonally ice-free solution (red curve), in contrast to the results 
reported by Thorndike. A second stable state which is perennially ice-free exists for 
the solution indicated by the black curve, as discussed by Thorndike, but not for the 
solution shown here by the blue curve. A second stable state which is perennially ice- 
free does however exist for the solution indicated by the red curve. The coexistence 
in Thorndike's model of a stable seasonally ice-free solution and a stable perennially 
ice-free solution is consistent with the results presented here (Fig. 3 of the Research 
Report). 
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